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LONG  TERM  GOALS 

My  ultimate  goal  in  this  work  is  to  develop  a  capability  and  physical  basis  to  simulate  oriented  fractures 
and  leads  in  sea  ice  in  a  way  that  can  be  utilized  in  the  next  generation  of  sea  ice  models. 

OBJECTIVES 

The  main  objective  of  this  effort  is  to  develop  and  numerically  investigate  a  physically  based  model  for 
describing  and  simulating  the  failure  and  flow  of  anisotropic  cracked  sea  ice,  and  ultimately  simulate 
oriented  fractures  on  the  geophysical  scale  in  sea  ice.  With  this  objective  in  mind  a  conceptual  model 
based  on  oriented  lead  ice  imbedded  in  thicker  ice  has  been  formulated  numerically  with  the  yield 
criterion  for  both  the  thick  and  thin  ice  taken  from  laboratory  observations  together  with  an  associated 
non  normal  flow  rule  prescribed  in  an  energetically  consistent  manner.  I  then  plan  to  utilize  this  model 
together  with  wind  and  water  stress  forcing  on  the  ice  pack  to  simulate  oriented  fractures  in  pack  ice 
and  compare  these  features  with  observations  from  satellite  imagery.  Secondary  related  objectives  are 
(1)  to  develop  a  capability  to  simulate  high  frequency  variability  in  pack  ice  as  this  should  prove  critical 
in  the  actual  formation  of  opening  leads  as  opposed  to  sliding  fractures  and  (2)  develop  an 
energetically  consistent  ice  dynamics  model  in  curvilinear  co-ordinates  to  be  used  both  in  forecast 
applications  and  for  consistently  determining  positive  definite  mechanical  energy  dissipation  for  use  in 
ambient  noise  acoustic  models  for  the  Arctic. 

APPROACH 

A  conceptual  notion  especially  relevant  to  the  geophysical  scale  is  that  leads  and  flaws  develop  from  the 
interaction  and  propagation  of  existing  oriented  flaws.  To  formulate  this  numerically,  the  fundamental 
model  consists  of  an  oriented  weak  lead/flaw  imbedded  in  thicker  ice  as  shown  in  Figure  1. 
Commensurate  with  continuum  mechanics  postulates  the  interfacial  stresses  (an  and  an)  along  the 
lead  are  taken  to  be  the  same  in  the  thick  and  thin  ice  while  stresses  perpendicular  to  the  lead  strike 
(022)  may  be  different.  For  the  basic  numerical  studies  the  isotropic  constitutive  laws  for  both  the  thick 
and  thin  ice  are  taken  to  be  coulombic  in  character  (Figure  1-b)  with  slopes  and  tensile  stresses 
consistent  with  laboratory  observations.  Other  constitutive  laws  are  however  also  utilized.  An 
associated  non  normal  flow  rule  for  this  yield  curve  is  constructed  to  insure  energetic  consistency  (i.e., 

the  quantity  Eus  jjajj  is  always  positive  definite).  In  addition  the  transition  from  coulombic  flow  to  out 
of  plane  failure  is  taken  to  occur  at  the  transition  between  diverging  and  converging  flow.  For 
numerical  investigation  this  model  is  generalized  to  include  oriented  leads  in  a  variety  of  directions.  A 
series  of  numerical  investigations  of  this  model  under  both  near  field  stress  forcing  and  far  field  stress 
are  carried  out.  In  the  far  field  stress  forcing  experiments  oriented  flaws  may  interact  with  other 
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oriented  flaws  at  different  locations.  These  results  are  analyzed  to  determine  basic  fracture  propagation 
characteristics  as  well  as  deformation  behavior  with  the  goal  of  simulating  oriented  fractures  and  leads. 


(a) 


Figure  1.  (a)  Schematic  of  oriented  thin  ice  lead  or  flaw  embedded  in  thick  ice.  An  isotropic 
elliptical  yield  curve  for  both  the  thin  and  thick  ice  is  shown  in  the  inset,  (b)  Laboratory  based 
coulombic  yield  curve.  Stresses  above  or  to  the  right  of  the  dotted  lines  are  tensile. 


WORK  COMPLETED 


A  series  of  numerical  experiments  using  several  basic  yield  curves  were  carried  out  in  order  to 
extensively  numerically  investigate  the  near  and  far  field  failure  characteristics  of  the  basic  model. 
Numerical  experiments  were  also  carried  out  to  determine  precisely  where  the  transition  between 
opening  and  closing  leads  occurs  and  the  necessary  conditions  for  energetic  consistency.  Both  analytic 
and  numerical  investigations  were  carried  out  to  elucidate  the  physical  principles  underlying  the 
modification  of  the  aggregate  coulombic  failure  for  this  model  even  though  there  is  a  constant  slope 
failure  in  the  basic  yield  curve.  To  understand  the  far  field  behavior  a  series  of  far  field  experiments 
were  carried  out  with  both  interacting  and  isolated  oriented  flaws  surrounded  by  thicker  ice  (Hibler  and 
Schulson,  1998).  To  develop  a  true  fracture  propagation  model  a  general  procedure  for  integrating  the 
yield  curve  was  initiated. 

With  respect  to  other  work,  a  procedure  for  imbedding  the  ice  in  the  boundary  layer  was  developed  for 
investigating  the  interaction  between  propagating  kinematic  waves  and  inertial  oscillations.  This  model 
(Hibler  et.  al.,  1998)  was  used  in  conjunction  with  observational  data  on  ice  deformation  to  determine  if 
this  mechanism  could  be  responsible  for  much  of  the  high  frequency  variability  that  is  a  ubiquitous 
feature  of  observed  ice  deformation. 


Finally,  a  differencing  procedure  for  curvilinear  co-ordinates  was  developed  that  insures  that  (a)  the 
globally  integrated  work  done  by  ice  interaction  forces  equals  the  global  mechanical  energy  dissipation, 
and  (b)  the  mechanical  energy  dissipation  is  always  positive  definite.  For  momentum  conservation  it 
was  found  that  this  differencing  was  considerably  different  than  that  used  in  rectangular  co-ordinates. 


RESULTS 


The  essential  result  is  that  with  an  anisotropic  model  of  a  single  flaw  preferred  lead  orientation  occurs 
under  near  field  stresses  for  ice  compressed  in  a  biaxial  state.  In  contrast  to  traditional  coulombic 
conjugate  plane  theory  of  faults  for  which  a  constant  orientation  would  be  expected,  the  preferred  lead 
orientation  depends  on  the  confinement  ratio  with  lower  confinement  ratios  producing  leads  more 
closely  aligned  with  the  largest  principal  stress  with  the  intersection  angle  going  to  zero  for  large 
enough  positive  confinements.  If  this  result  is  formulated  in  terms  of  an  effective  friction  coefficient  the 
results  show  a  decrease  in  the  effective  friction  for  larger  confinement,  a  result  consistent  with 
observations  in  rock  mechanics.  Analysis  of  the  results  indicates  that  these  features  arise  because  of  the 
finite  width  of  the  lead  allowing  compressive  stress  along  the  lead. 

While  for  a  given  confinement  ratio  the  precise  orientation  of  the  leads  does  depend  on  the  constitutive 
law  assumed  for  the  thin  and  thick  ice,  if  we  consider  a  flaw  imbedded  in  surrounding  stronger  ice  and 
apply  a  far  field  stress  to  that  region  the  behavior  is  quite  different.  In  this  instance  we  find  that  the  local 
stress  field  is  modified  by  the  presence  of  the  oriented  flaw  which  in  turn  modifies  the  nearby  maximum 
tensile  stresses  (Figure  2).  In  the  case  of  different  flaw  orientations  we  find  that  there  is  a  much 
narrower  range  of  orientations  (relative  to  the  far  field  stress)  that  yield  the  greatest  tensile  stresses  in 
the  neighborhood  of  the  flaw.  In  a  manner  relatively  independent  of  strength,  flaws  oriented  at  about 
20*  relative  to  the  far  field  stresses  tend  to  yield  the  most  pronounced  tensile  stresses.  As  shown  in 
Figure  2  this  characteristic  is  somewhat,  but  not  significantly, 

affected  by  the  isotropic  rheology  assumed  for  the  thin  and  thick  ice.  These  results  and  theoretical 
model  are  described  in  a  journal  article  by  Hibler  and  Schulson  (1998). 


Maximum  Tensile  Stress 
Adjacent  to  Flaw 


Figure  2.  (a)  Maximum  tensile  stress  adjacent  to  flaw  as  a  function  of  angle  of  the  flaw  relative  to 
the  far  field  forcing.  Curves  for  different  strengths  of  flaws  are  shown  with  an  elliptical  yield  curve 
as  shown  in  the  inset  of  Figure  la.  The  laboratory  yield  curve  is  a  mohr  coulombic  type  yield  curve 
as  shown  in  Figure  lb.  (see  Hibler  and  Schulson,  1998).  (b)  Orientation  of  embedded  flaw 
surrounded  by  isotropic  ice  under  far  field  forcing. 


These  far  field  results  suggests  a  different  mechanism  than  conjugate  coulomb  faulting  from  a  parent 
isotropic  coulomb  rheology  is  responsible  for  lead  formation.  With  this  far  field  result  in  mind  we  can, 
for  example,  imagine  that  a  complex  stress  field  existing  in  pack  ice  (due  predominately  to  wind  and 
water  stresses)  will  tend  to  favor  the  weakening  and  further  propagation  of  flaws  with  orientations 
roughly  at  20*  relative  to  the  'far  field'  pack  ice  stress.  Other  leads  will  tend  to  close  and  strengthen  or 
weaken  more  slowly.  Because  of  this  selection  phenomenon  a  good  first  approximation  to  this 
phenomenon  may  be  gotten  by  beginning  with  a  uniform  array  of  oriented  flaws  in  every  grid  cell.  The 
complex  stress  field  will  eventually  lead  to  the  failure  of  two  intersecting  leads  at  some  point,  only  one 
of  which  will  survive  a  loading  of  the  stress  field  because  of  the  unstable  character  of  the  two  flaw 
system.  Once  any  oriented  fracture  forms  in  this  manner  the  system  then  becomes  locally  anisotropic. 
As  the  stress  state  increases  or  the  formed  flaw  weakens  by  opening,  the  nearby  tensile  stresses  will 
activate  other  adjacent  pairs  of  oriented  leads.  As  before,  only  one  of  these  leads  will  survive  and  the 
fracturing  process  will  continue  with  the  failure  not  only  occurring  over  a  narrow  region  but  also 
containing  oriented  flaws. 


To  examine  high  frequency  variability  issues  ,  an  idealized  1.5  dimensional  dynamic-thermodynamic  sea 
ice  model  was  configured  to  be  relevant  to  results  of  a  winter  deformation  experiment  in  a  tide  free 
region  of  the  East  Antarctic  sea  ice  zone  (Worby  and  others,  1996).  During  this  time,  a  number  of 
drifting  buoys  were  monitored  hourly  for  studies  of  sea  ice  drift  and  deformation  (Heil  and  others, 

1998)  with  results  yielding  substantial  inertial  signals  in  the  deformation.  The  model  configuration 
consisted  of  a  40km  lagrangian  grid  with  wind  forcing  taken  to  be  perpendicular  to  the  coast.  In  this 
model  (Hibler  et.  al.,  1998),  full  two  dimensional  ice  dynamics  is  assumed  but  all  variations  along  the 
coast  are  assumed  to  be  zero.  For  simplicity  pressure  gradient  variations  in  the  ocean  due  to  Ekman 
convergence  were  also  neglected. 

Simulations  with  this  grid  with  different  wind  forcing  show  that  without  an  improved  boundary  layer 
formulation  both  the  ice  drift  and  deformation  follow  the  wind  forcing  yielding  no  inertial  deformation. 
With  an  improved  boundary  layer  formulation  but  no  ice  thickness  evolution  equations  ice  drift  does 
show  an  inertial  signal,  but  there  is  no  inertial  deformation  signal,  except  at  one  grid  cell  right  at  the 
coast.  With  thickness  strength  coupling,  however,  there  is  a  substantial  signal  in  the  ice  deformation 
with  kinematic  waves  (see  e.g.,  Hibler  et.  al.,  1983)  propagating  outward  interacting  with  the  inertial 
motions  to  produce  varying  deformation  in  regions  far  from  the  coast. 

These  oscillations  are  notably  present  even  in  the  absence  of  wind  as  shown  in  Figure  3.  In  this  figure 
we  show  deformation  results  after  a  spatially  constant  pulse  wind  was  applied  to  the  system  for  about  4 
hours  and  then  dropped  to  zero.  As  can  be  seen  there  is  a  substantial  change  in  the  ice  concentration 
with  eventually  a  "banded"  structure  developing  with  oscillating  regions  of  thin  ice  alternating  spatially 
with  ice  undergoing  no  deformation.  Because  of  the  non-linear  sea  ice  dynamics  there  is  only  a  weak 
damping  of  these  oscillations  although  they  will  eventually  die  out  leaving  "banded"  ice  concentrations 
behind.  Comparison  of  the  magnitude  of  these  oscillations  with  those  observed  during  quiescent  wind 
conditions  (Heil  et.  al.,  1998)  show  them  to  be  very  close  in  magnitude. 

ICE  COHCBnRmOM  UC  IK  HD  WCTHtCE  HUM  CHORE 


Figure  3.  Perspective  view  of  compactness  as  a  function  of  time  and  distance  from  shore 
about  fifty  hours  after  an  initial  pulse  of  constant  onshore  wind  stress. 


Finally  with  regard  to  mechanical  energy  dissipation  in  sea  ice  dynamics  models,  an  energetically 
consistent  finite  difference  formulation  in  spherical  co-ordinates  was  developed.  Tests  with  the 
cartesian  co-ordinate  version  of  this  model  (Hibler  and  Ip,  1995)  showed  that  the  addition  of  energetic 
consistency  led  to  a  substantial  improvement  in  the  agreement  between  simulated  and  observed  buoy 
drift,  especially  as  regards  excessive  stoppage  of  ice  motion.  Formulation  of  the  energetic  consistent 
equations  in  curvilinear  co-ordinates  was  found  to  require  that  all  metric  terms  be  incorporated  into  the 
differencing  scheme.  Analytical  proofs  of  the  differencing  scheme  were  constructed  which  showed 
explicitly  both  momentum  and  energy  conservation.  Numerical  results  together  with  documentation  of 
the  model  and  development  of  the  continuum  equations  are  included  in  a  report  (Hibler,  1998)  prepared 
under  the  auspices  of  the  Arctic  Chair  and  the  Naval  Postgraduate  School  for  inclusion  in  ongoing  PIPS 
modeling  efforts. 

IMPACT/APPLICATIONS 

These  results  are  beginning  to  have  an  impact  on  ice  forecasting  especially  within  the  context  of 
mechanical  energy  output  for  use  in  acoustic  applications.  The  theoretical  development  of  anisotropic 
models  are  also  impacting  both  laboratory  and  field  observations  and  thoughts  on  the  basic  mechanisms 
of  lead  and  fracture  formation. 

TRANSITIONS 

The  main  transition  is  the  development  and  documentation  of  a  fully  energy  conserving  curvilinear  co¬ 
ordinate  model  for  use  in  PIPS  3.0  and  in  ambient  noise  acoustic  applications. 

RELATED  PROJECTS 

Closely  related  to  this  project  is  one  by  Erland  Schulson  under  a  separate  ONR  grant  N00014-97-1- 
0211  on  laboratory  experiments  on  both  the  small  and  intermediate  scale  to  observe  mechanisms  and 
characteristics  of  ice  fracturing  and  failure.  The  synergy  between  these  efforts  on  the  mechanisms  for 
fracture  have  substantially  aided  in  the  modeling  work  and  selection  of  numerical  experiments. 

REFERENCES 

Heil,  P.,  V.  I.  Lytle  and  I.  Allison,  1998:  Enhanced  thermodynamic  ice  growth  by  sea  ice  deformation, 
Ann.  Glaciol.,  27,  in  press. 

Hibler,  W.  D.  Ill,  and  E.  M.  Schulson,  1998:  On  a  Model  for  the  Anisotropic  Failure  of  Flawed  Sea 
Ice,  Journal  of  Geophysical  Research,  accepted  for  publication. 

Hibler,  W.  D.  Ill,  P.  Heil  and  V.  Lytle,  1998:  Modeling  High  Frequency  Variability  in  Antarctic  Sea  Ice 
Models,  Annals  of  Glaciology,  27,  in  press. 

Hibler,  W.  D.,  Ill,  I.  Udin,  and  A.  Ullerstig  (1983),  On  forecasting  mesoscale  ice  dynamics  and  buildup, 
Annals  of  Glaciology,  4,  110-115. 


Hibler,  W.  D.  Ill,  and  C.  F.  Ip,  1995:  The  effect  of  sea  ice  rheology  on  Arctic  buoy  drift,  in  Ice 
Mechanics,  ASME  AMD-Vol.  207,  Dempsey  and  Rajapakse,  Eds.,  Book  No.  H00954-1995,  255-263. 


Hibler,  W.  D.  Ill,  1998:  A  curvilinear  co-ordinate  energy  conserving  sea-ice  dynamics  model  for  PIPS, 
NPS  report,  in  preparation. 

Worby,  A.  P.,  N.  L.  Bindoff,  V.  I.  Lytle,  I.  Allison  and  R.  A.  Massom,  1996:  Winter  ocean/sea  ice 
interactions  studies  in  the  East  Antarctic,  EOS  Transactions  77(46),  456-457. 

PUBLICATIONS 

Hibler,  W.  D.  Ill  and  E.  M.  Schulson,  1997:  On  modeling  sea-ice  fracture  and  flow  in  numerical 
investigations  of  climate,  Annals  of  Glaciology,  25,  26-32,  1997. 

Hibler,  W.  D.  Ill,  and  E.  M.  Schulson,  1998:  On  a  Model  for  the  Anisotropic  Failure  of  Flawed  Sea 
Ice,  Journal  of  Geophysical  Research,  accepted  for  publication. 

Hibler,  W.  D.  Ill,  P.  Heil  and  V.  Lytle,  1998:  Modeling  High  Frequency  Variability  in  Antarctic  Sea  Ice 
Models,  Annals  of  Glaciology,  27,  in  press. 


